import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import kstest

# 设置字体和显示格式，确保中文显示正确
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用SimHei字体
plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号

# 读取Excel文件，保留原始列名
filename = '附件2_2.xlsx'  # 替换为实际文件路径
data = pd.read_excel(filename)

# 定义函数来处理字符串中的范围并取其均值
def convert_to_mean(value):
    if isinstance(value, str) and '-' in value:
        parts = value.split('-')
        try:
            # 将范围内的值取均值
            numbers = [float(part) for part in parts]
            return np.mean(numbers)
        except ValueError:
            return np.nan  # 如果转换失败则返回 NaN
    else:
        try:
            return float(value)
        except ValueError:
            return np.nan  # 如果不是数值则返回 NaN

# 对需要转换的数据列进行处理，确保数据为数值形式
data['销售单价/(元/斤)'] = data['销售单价/(元/斤)'].apply(convert_to_mean)

# 提取数据并删除缺失值
yield_per_acre = data['亩产量/斤'].dropna()  # 删除亩产量中的缺失值
cost_per_acre = data['种植成本/(元/亩)'].dropna()  # 删除种植成本中的缺失值
price_per_kg = data['销售单价/(元/斤)'].dropna()  # 删除销售单价中的缺失值

# 定义函数：正态性检验，使用K-S检验
def is_normal_distribution(data):
    # 如果P值大于0.05，则认为数据符合正态分布
    return kstest((data - np.mean(data)) / np.std(data), 'norm')[1] > 0.05

# 检查三个变量是否符合正态分布
is_normal_yield = is_normal_distribution(yield_per_acre)
is_normal_cost = is_normal_distribution(cost_per_acre)
is_normal_price = is_normal_distribution(price_per_kg)

# 可视化：创建一个3x1的子图，用于展示异常值处理情况
fig, axs = plt.subplots(3, 1, figsize=(10, 12))

# 绘制亩产量的异常值处理图
if is_normal_yield:
    # 如果亩产量符合正态分布，使用3σ原则检测异常值
    mu_yield = np.mean(yield_per_acre)
    sigma_yield = np.std(yield_per_acre)
    outliers_yield = yield_per_acre[np.abs(yield_per_acre - mu_yield) > 3 * sigma_yield]

    # 绘制直方图并标出3σ区域和异常值
    axs[0].hist(yield_per_acre, bins=30, density=True, alpha=0.6, color='g')
    axs[0].axvline(mu_yield - 3 * sigma_yield, color='r', linestyle='--', label='3σ下限')
    axs[0].axvline(mu_yield + 3 * sigma_yield, color='r', linestyle='--', label='3σ上限')
    axs[0].scatter(outliers_yield.index, outliers_yield, color='r', label='异常值')
    axs[0].set_title('亩产量 (正态分布 - 3σ原则)')
else:
    # 如果亩产量不符合正态分布，使用箱型图来显示异常值
    axs[0].boxplot(yield_per_acre)
    axs[0].set_title('亩产量 (非正态分布 - 箱型图)')

# 绘制种植成本的异常值处理图
if is_normal_cost:
    mu_cost = np.mean(cost_per_acre)
    sigma_cost = np.std(cost_per_acre)
    outliers_cost = cost_per_acre[np.abs(cost_per_acre - mu_cost) > 3 * sigma_cost]

    axs[1].hist(cost_per_acre, bins=30, density=True, alpha=0.6, color='g')
    axs[1].axvline(mu_cost - 3 * sigma_cost, color='r', linestyle='--', label='3σ下限')
    axs[1].axvline(mu_cost + 3 * sigma_cost, color='r', linestyle='--', label='3σ上限')
    axs[1].scatter(outliers_cost.index, outliers_cost, color='r', label='异常值')
    axs[1].set_title('种植成本 (正态分布 - 3σ原则)')
else:
    axs[1].boxplot(cost_per_acre)
    axs[1].set_title('种植成本 (非正态分布 - 箱型图)')

# 绘制销售单价的异常值处理图
if is_normal_price:
    mu_price = np.mean(price_per_kg)
    sigma_price = np.std(price_per_kg)
    outliers_price = price_per_kg[np.abs(price_per_kg - mu_price) > 3 * sigma_price]

    axs[2].hist(price_per_kg, bins=30, density=True, alpha=0.6, color='g')
    axs[2].axvline(mu_price - 3 * sigma_price, color='r', linestyle='--', label='3σ下限')
    axs[2].axvline(mu_price + 3 * sigma_price, color='r', linestyle='--', label='3σ上限')
    axs[2].scatter(outliers_price.index, outliers_price, color='r', label='异常值')
    axs[2].set_title('销售单价 (正态分布 - 3σ原则)')
else:
    axs[2].boxplot(price_per_kg)
    axs[2].set_title('销售单价 (非正态分布 - 箱型图)')

# 添加全局标题
plt.suptitle('作物亩产量、种植成本和销售单价的异常值检测')

# 调整布局并显示图表
plt.tight_layout(rect=[0, 0, 1, 0.95])  # 调整图表和标题的间距
plt.savefig('fig/异常值检测结果.png')  # 保存图表为PNG文件
plt.show()

# 绘制QQ图
from scipy import stats
fig, axs = plt.subplots(3, 1, figsize=(10, 15))

# 亩产量QQ图
stats.probplot(yield_per_acre, dist="norm", plot=axs[0])
axs[0].set_title('亩产量 QQ图')

# 种植成本QQ图
stats.probplot(cost_per_acre, dist="norm", plot=axs[1])
axs[1].set_title('种植成本 QQ图')

# 销售单价QQ图
stats.probplot(price_per_kg, dist="norm", plot=axs[2])
axs[2].set_title('销售单价 QQ图')

plt.tight_layout()
plt.savefig('fig/qq_plots.png')  # 保存图表为PNG文件
plt.show()
